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We develop analytic approximations of thermodynamic functions of fully ionized nonideal electron- 
ion plasma mixtures. In the regime of strong Coulomb coupling, we use our previously developed 
analytic approximations for the free energy of one-component plasmas with rigid and polarizable 
electron background and apply the linear mixing rule (LMR). Other thermodynamic functions are 
obtained through analytic derivation of this free energy. In order to obtain an analytic approximation 
for the intermediate coupling and transition to the Debye-Huckel limit, we perform hypernetted- 
chain calculations of the free energy, internal energy, and pressure for mixtures of different ion 
species and introduce a correction to the LMR, which allows a smooth transition from strong to 
weak Coulomb coupling in agreement with the numerical results. 

PACS numbers: 52.25.Kn, 05.70.Ce, 52.27.Gr 



I. INTRODUCTION 



II. PLASMA PARAMETERS 



We study the equation of state (EOS) of fully ionized 
nonideal electron-ion plasmas (EIP). In a previous work 
[E 01) hypernetted chain (HNC) calculations were per- 
formed and analytic formulae were proposed for EOS 
calculations of EIP containing a single ion species. For 
mixtures of different ion species, the EOS was calculated 
using the linear mixing rule (LMR) , whose high accuracy 
at r > 1 was previously confirmed in a number of studies 
[E 0) IE IE • However, the LMR is inaccurate for weakly 
coupled plasmas. Some consequences of its violation were 
studied by Nadyozhin and Yudin Q, who showed that 
the differences between the linear and nonlinear mixing 
at moderate Coulomb coupling (0.1 < Y < 1) can shift 
the nuclear statistical equilibrium at the final stage of a 
stellar gravitational collapse. 



In this paper, we perform HNC calculations of the free 
energy, internal energy, and pressure for mixtures of var- 
ious kinds of ions in the weak, intermediate, and strong 
coupling regimes, and suggest an analytic correction to 
the LMR. 



In Sec. |TT] we define the basic plasma parameters. In 
Sec. IIIII we calculate the EOS of ion mixtures and propose 
an analytic formula for the EOS of multicomponcnt EIP, 
applicable at any T values. The summary is given in 
Sec. HV1 



Let n e be the electron number density and rij the num- 
ber density of ion species j=l,2,. . . , with mass and charge 
numbers Aj and Zj , respectively. The total number den- 
sity of ions is rii on = Y], rij . The electric neutrality 
implies n e — (Z)n[ on . Here and hereafter the brackets 
(. . .) denote averaging of the type (/) = Xjfj, where 
Xj — Ti j /n[ on . 

The state of a free electron gas is determined by the 
electron number density n e and temperature T. Instead 
of n e it is convenient to introduce the dimensionless den- 
sity parameter r s — a e /ao, where a e = (|7m e ) -1 / 3 and 
ao is the Bohr radius. 

At stellar densities it is convenient to use, instead of r s , 
the relativity parameter [l(| x rc \ = pp/m e c = 0.014 rj 1 , 
where pf — h (37r 2 rt e ) 1 / 3 is the electron Fermi momen- 
tum. The Fermi kinetic energy is €f — c \J (m e c) 2 + p F — 
m e c 2 , and the Fermi temperature equals Tf = ep/ks = 
T r {lr - 1), where T r = m e c 2 /k B = 5.93 x 10 9 K, 
7 r = yT + x' 2 cl , and ks is the Boltzmann constant. 

The ions are nonrelativistic in most applications. The 
strength of the Coulomb interaction of ion species j is 
characterized by the Coulomb coupling parameter, 



= (Z ie ) 2 /( aj k B T) = T e Z 

A/3 



5/3 



(1) 
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where aj = a e Zj is the ion sphere radius and T e = 
e 2 j (a e kBT). In a multicomponent plasma, it is useful to 
introduce the mean ion-coupling parameter T = T e (Z 5 ^ 3 ) 
3]. At a melting temperature T m , corresponding toTsj 
175 (e.g., Q), the plasma freezes into a Coulomb crystal. 

An important scale length is the thermal de Broglie 
wavelength \j — (2t:H 2 /m-jkBT) 1 / 2 , where rrij is the ion 
mass. The electron thermal length A e is given by the 
same expression with rrij replaced by m e . 

The quantum effects on ion motion become impor- 
tant at T < T p , where T p = hajp/kg and lu p = 



2 



(47re 2 nj on (Z 2 /m)) is the ion plasma frequency. 

In this paper we consider only the classical Coulomb 
liquid, which implies T ^> T p and T > T rn . 



analytic approximation [![. The electron-gas contribu- 
tions to xt, CV> and S tend to zero at T <C Tp. Related 
numerical problems and their cure will be discussed else- 
where (l2| . 



III. EQUATION OF STATE 



Assuming commutativity of the kinetic and potential 
operators and separation of the traces of the electronic 
and ionic parts of the Hamiltonian, the total Helmholtz 
free energy F can be conveniently written as 



F = FT 



where F™ and F^ denote the ideal free energy of ions 
and electrons, and the last three terms represent an ex- 
cess free energy arising from the electron-electron, ion- 
ion, and ion-electron interactions, respectively. 

The pressure P, the internal energy U, and the entropy 
S of an ensemble of fixed number of plasma particles 
in volume V can be obtained using the thermodynamic 
relations P = -(dF/8V) T , S = -(dF/dT) v , and 
U = F+TS. The second-order thermodynamic functions 
are derived by differentiating these first-order ones. The 
decomposition j5J induces the analogous decomposition 
of P, U, S, the heat capacity Cy — (dS/dlnT)y, and the 
logarithmic pressure derivatives xt — (<91nP/<91nT)y 
and Xp = —{d\nP/d\nV)T- Other second-order func- 
tions can be expressed through these by Maxwell rela- 
tions (e.g., (TT1|). 



F 



(e) 



id 



F R , 



Fa ~\~ Fit 



(2) 



A. Ideal electron-ion plasmas 

The free energy of a gas of Nj — njV nonrelativistic 
classical ions of jth kind is 



F$ = NjksT [infaAf/fc) 



1] 



(3) 



where gj is the spin multiplicity. The total free energy 
is given by the sum F; 1 ™ = J2j ^la ■ Analogous sums 



B. Nonideal plasmas containing one type of ions 

Let us recall fit formulae for nonideal EIP containing 
a single kind of ions. 

a. Electron exchange and correlation. Electron- 
electron (exchange-correlation) effects were studie d by 
many authors. For the reasons explained in Refs. [llll2l|. 
we adopt the fit to f ee = F ee /(N e kBT) presented in 
Ref. pj|. 

b. One- component plasma. The internal energy of 
the liquid one-component plasma (OCP) at any values 
of r is given by [|J 



= r 3/2 



At 



-4, 



.VM+V i + r 

where Uu = Uu/kBTN- wu , and 



gi r 2 



A* 



-V3/2 



p2 ' 

(5) 



(6) 



ensures the correct transition to the Debye-Hiickel limit. 
The parameters A x = -0.907347, A 2 = 0.62849, B x = 
0.0045, B 2 = 170, B 3 = -8.4 x 10~ 5 , and B 4 = 0.0037 
allow one to reproduce the best available MC simulations 
of liquid OCP at 1 < V < 190 [3 with an accuracy 
matching the numerical MC noise. 

From Eq. ([5]) one obtains the analytic expression for 
fa = Fu/kBTNi on by integration, and then the Coulomb 
contributions to the other thermodynamic functions by 
differentiation (2j. 

c. Electron polarization. Electron polarization in 
Coulomb liquid was studied by perturbation 0, 
and HNC 0, 0, HJ techniques. The results for_/ i( 



o la t-^j ia o e— > 1 — u r—L" 

give U, S, P, and C v . Since Eq. © contains rij under F le /N ion k B T have been fitted by the expression [2| 



logarithm, these sums for F and S naturally include the 
entropy of mixing S'mix = —ks Ylj m a;j. 
The free energy of the electron gas is given by 



= ^N e 



P-\ e > V 

1 iH ' 



(4) 



where /i e is the electron chemical potential. The pressure 
P^ and the number density n e = N e /V are functions of 
He and T, which can be written through the Fcrmi-Dirac 
integrals I v (xei T )i where %e = He/ksT and v = 1/2, 
3/2, and 5/2. In Ref. [l| we gave analytic approxima- 
tions for the Fermi-Dirac integrals, based on fits [9j] to 
electron-positron thermodynamic functions. The chemi- 
cal potential at a given density can be found either from 
a numerical inversion of function n e {x ei T) or using the 



, _ _p c DH /rI+ c TF ar^g 1 (r s ,T e )g 3 (x TC i) 

6 1+ [6^ + a 52 (r s ,r e )r-/r s ] 7l 7 1, U 

The coefficients cdh, ctf, a, b, v and functions 
51,2(^8) Te) and g^(x Te \) parametrically depend on the ion 
charge Z . Here the coefficients Cdh and ctf are not free 
fit parameters, because 



c DH = (Z/y/3) [(1 + Zf' 2 - 1 - Z 3 / 2 } 



(8) 



ensures the transition of the excess free energy to the 
Debye-Hiickel limit at small T, and ctf at large Z is 
given by the Thomas- Fermi theory (Toj . 




FIG. 1: Fractional difference between the Coulomb part of 
the internal energy (Ua) in different approximations and the 
LMR prediction as a function of the average ion Coulomb cou- 
pling parameter F for binary mixtures of ions with Z2 fZ\ = 2 
and X2 — 1 — x\ = 0.05 and 0.2 in the rigid background. 
Dashed lines (DH): Debye-Hiickel formula, dot-dashed lines 
(CLM): corrected linear mixing dots: MC results of 
Ref. [y] for X2 = 0.05; crosses: HNC results of Ref. trian- 
gles: present HNC results; solid lines: present fit. 



C. Nonideal mixtures of ions 

A common approximation for the excess (nonideal) free 
energy of the strongly coupled ion mixture is the LMR, 



fax (r) ~ / ' x jfcx(^j, Xj — 1) . 



(9) 



where superscript "LM" denotes the linear-mixing ap- 
proximation, and all Tj correspond to the same T e (as- 
suming that the pressure is given almost totally by the 
strongly degenerate electrons): Tj = T Z^ 3 / (Z 5 / 3 ) . In 
Eq. ([9|), /ox is the reduced nonideal part of the free 
energy: / cx = fa for the "rigid" (uniform) charge- 
neutralizing electron background and f cx = fa + he + 
Zjfee for the polarizable background 

The high accuracy of Eq. ([9]) for binary ionic mixtures 
in the rigid background was first demonstrated by cal- 
culations in the HNC approximation 0] and confirmed 
later by MC simulations (e.g., 0, @, 0] ) ■ 

The validity of the LMR in the case of an ionic mix- 
ture immersed in a polarizable finite-temperature elec- 
tron background has been examined by Hansen et al. Q 
in the first-order thermodynamic perturbation approxi- 
mation and by Chabrier and Ashcroft Q by solving the 
HNC equations with effective screened potentials. These 



FIG. 2: Fractional differences between the nonideal part of 
the free energy due to the ion-ion and ion-electron interactions 
(Fu+Fie) in different approximations and the LMR prediction 
for a mixture of ions with Zi = 1 and Z2 = 8 and X2 = 
1 — xi — 0.25 in the rigid and polarizable electron background, 
as functions of F. Dashed lines: DH, dot-dashed lines: CLM 
asterisks: HNC results of Ref. Q; solid lines: present fit. 



authors found that the LMR remains accurate when the 
electron response is taken into account in the inter-ionic 
potential, as long as the Coulomb coupling is strong 
(T > 1). 

On the other hand, the LMR is invalid at F <C 1. 
Indeed, in this case the Debye-Hiickel theory gives 



Mi ' 



(10) 



fDH _ _p3/2 1/0 fDH _ f DH/I 
Jee — 1 e / v °i J it ~ Jee Sil 

and for the EIP 

fDH _ fDH 1 fDH 1 / 7 \ f DH_ /-DHaDH /-, -, \ 

Jcx - hi + He + \ /j lJee — Jee Scip 1 l 1J J 

where 

,DH = (Z 2 ) 3/2 ,DH = ((Z 2 ) + (Z)) 3 ^ 

However, the LMR at F <C 1 gives another result: 

fLM ^ fDH/LM fLM ^ fDH aLM /-i o\ 

J ii Jee ^ii ' J ex Jee ^eip ' V / 



where 



C^ = (Z^ 2 ), C^ = (Z(Z + l)^). (14) 



Nadyozhin and Yudin [8| considered several possible 
modifications of the LMR at intermediate T and found 
that such modification can appreciably shift the statis- 
tical nuclear equilibrium at the conditions typical of the 
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FIG. 3: Fractional difference between the Coulomb part of 
the internal energy in different approximations and the LMR 
prediction for a binary mixture of ions with Z2 jZ\ = 8 in the 
rigid background. Dashed lines (DH): Debye-Hiickel formula, 
symbols: present HNC results; solid lines: present fit. Differ- 
ent symbols correspond to different xi values: 0.01 (squares), 
0.05 (solid triangles), 0.2 (dots), 0.25 (asterisks), 0.5 (empty 
triangles), and 0.7 (empty circles). 



FIG. 4: Fractional difference between the Coulomb part of 
the free energy and the LMR prediction for binary ionic mix- 
tures with ZijZ\ = 12 and Z2/Z1 = 16, and with X2 = 0.01, 
0.05, and 0.2. Dashed lines (DH): Debye-Hiickel formula, 
symbols: present HNC results; solid lines: present fit for 
Z2/Z1 = 16; long-dash-dot lines: present fit for Z2/Z1 — 12. 
Different symbols correspond to different combinations of X2 
and Z2/Z-1 values. 



final stage of a stellar gravitational collapse. They con- 
sidered the rigid background and advocated a modifica- 
tion of every term in Eq. @ by multiplying the leading 
fit coefficient at small T by factor dj = \J (Z 2 ) /Zj(Z). 
It corresponds to replacing V3/2 in Eq. ^ by dj y/3/2 
(and a possible simultaneous change of A%). For a com- 
pressible background an analogous modification implies 
additionally replacement of Eq. ([5]) by 



JJ) _ Z 3 

\/3 



3/2 



(Z) 3 / 2 



- 1 



(15) 



Hereafter modifications of this type will be called cor- 
rected linear mixing (CLM). 

The result of such modifications is shown by dot- 
dashed lines in Figs. Q] and [21 where we plot the ra- 
tios Amh/u™ and A/ ex // e L x M as functions of T. Here 
A/ = / - / LM and Au = u - u LM are the deviations of 
the reduced free and internal energies, respectively, from 
the LMR. In Fig. Q] the electron response is neglected 
(rigid background). We see that, for example, at T 10 
the CLM prescription gives Au and A/ corrections of 
about 1% in Fig.Q]and several percents in Fig. [2 whereas 
they must be much smaller according to the MC results 
[sill- Moreover, Au and A/ in the CLM approximation 
have the incorrect sign at T > 1 (note that according to 



Ref. [|| Au/u and A/// are negative at any T). Addi- 
tional modifications of the coefficient A2 in Eq. ^ also 
do not solve the problem. 

An alternative to the CLM, named "complex mixing" 
in Ref. Q, maintains the correct sign of Af and Au, but 
leads to still larger absolute values of these corrections 
(i.e., still slower recovery of the LMR) at large T. 

In order to find a more accurate approximation, we 
have performed HNC calculations of Fa, Uu and Pa for 
binary ionic mixtures in the rigid background for a broad 
r range and various values of the charge number ratios 
Z2/Z1 and fractional abundances X2 — 1 — X\. Some of 
the results are shown by triangles in Figs. [1] and O In 
agreement with the previous studies (e.g., [H), our nu- 
merical results show monotonically decreasing fractional 
deviations from the LMR with increasing T. The re- 
sults agree to at least 5 digits with those published in 
@ (crosses in Fig. [T]). At V > 1, HNC results tend to a 
constant residual within 1%, which is due to the intrinsic 
inaccuracy of the HNC approximation for strongly cou- 
pled plasmas because of the lack of the bridge functions 
in the diagrammatic representation of this approxima- 
tion. To prove this statement, in Fig. [Tj we plot by dots 
the values of Aita/uf^ 1 from MC simulations 0]. The 
latter simulations give tiny deviations from the LMR at 
r > 3, which are invisible in the figure scale. 
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FIG. 5: Fractional difference between the Coulomb part of 
the internal energy and the LMR prediction for binary ionic 
mixtures at Fi = 0.03 as function of the ion charge ratio 
Z%/Z\. Dashed lines (DH): Debye-Hiickel formula, symbols: 
present HNC results; solid lines: present fit. Different symbols 
correspond to different %i values: X2 — 0.01 (triangles), 0.05 
(dots), and 0.2 (asterisks). 



Asterisks in Fig. [5] correspond to the HNC calculations 
for mixtures in the polarizable electron background 
which give qualitatively the same results as the calcula- 
tions for the rigid background. The second dot-dashed 
curve in this figure shows CLM for the polarizable back- 
ground, according to Eq. (fT5|) . 

A correction to the linear mixing rule, which exactly 
recovers the Debye limit at T — > and the LMR at T ^> 1 , 
and which agrees with the HNC data, can be expressed 
by the following analytic fitting formula: 



A/ = / ( 



DH 



(1 + aT b ) c ' 



(16) 



where £ = Cu or Ceip for the rigid or polarizable back- 
ground, respectively, and the parameters a, b, and c de- 
pend on plasma composition as follows: 

2.2S+175 4 C LM ~ C DH 

a ~ 1-b ' (Z 5 / 2 ) ' 

b = d-° 2 , c=l + rf/6, d=(Z 2 )/(Z) 2 . 

Thanks to the simple form of this formula, its deriva- 
tives are also rather simple. For example, the corrections 
to the reduced internal energy and heat capacity read, 



FIG. 6: Fractional difference between the Coulomb part of 
the free energy and the LMR prediction for mixtures of three 
(solid lines) and two (dot-dashed and dotted lines) types of 
ions with charge ratios Zj/Zi = 2, 8, and 12. Different sym- 
bols correspond to HNC results with different combinations 
of Xj and Zj. Solid lines: Z2/Z1 = 8 or 12 (Z2 is marked 
near the curves, assuming Z\ = 1), X2 = 0.2; Z3/Z1 = 12, 
xs = 0.01. Dot-dashed lines: the same but without the third 
kind of ions (2:3 = 0). Dotted line: binary mixture with charge 
ratio Z3/Z1 — 2 (x2 = 0) and number fraction X3 = 0.2. 



respectively, 



Au 



Ac = 



abcT" 
! ~ 1 + aT b 
abcY b 



A/, 



(17) 



i + ar fc 



1W/ &2CI ? A/. (18) 
2 J (l + aL b ) 2 J V ' 



This approximation has been compared with our HNC 
calculations for binary ionic mixtures in the rigid elec- 
tron background at Z^/Z\ = 2, 5, 8, 12, and 16, with 
X2 = 1 — xi ranging from 0.01 to 0.7. Some of the results 
are shown in the figures. In all figures the Debye-Hiickel 
approximation is drawn by dashed lines for comparison. 
In Figs. Q] and discussed above, and in Figs. [3l [5] our 
approximation is drawn by solid lines. Figure [3] shows 
the fractional difference between the different results for 
internal energy and the LMR approximation for binary 
ionic mixtures with ZijZ\ = 8 and X2 = 0.01, 0.05, 0.2, 
0.25, 0.5, and 0.7 (marked near the curves). In Fig. [4] 
we show the case of the highest considered charge asym- 
metries: long-dash-dot lines show the approximation Q16p 
for Z2/Z1 = 12 and solid lines for Z^jZ\ = 16. Finally, in 
Fig. [5] we show a dependence of A fa / on the charge 
ratio Z 2 /Zi at Li = 0.03. Here different symbols corre- 
spond to the HNC results for different values of X2- 
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FIG. 7: Same as in Fig. [6] but for different charge ratios and 
fractional abundances. Solid lines: Z2/Z1 = 2 and 5 (marked 
near the curves), 22 = 0.05; Zz/Z\ — 12, xz = 0.01. Dot- 
dashed lines: the same but with xz — 0. Dotted line: the 
same with X2 = (2:3 = 0.01). 



We have also performed calculations for mixtures of 
ions of three different types on the rigid electron back- 
ground and compared the results with Eq. (flB|) . The re- 
sults of the comparison are shown in Figs. [5] and [7] Solid 
lines show the difference of the Coulomb free energy from 
the LMR according to Eq. ([TB)1 for 3-component mix- 
tures; for comparison, dot-dashed and dotted lines are 
plotted for 2-component mixtures with the same charge 
ratios; symbols represent the HNC results. In all consid- 
ered cases, adding a third component to a binary mixture 
increases the deviations of Eq. p6[) from HNC results by 
less than a factor of 1.5. We conclude that the agreement 
between the fit and numerical results remains satisfac- 
tory. 

IV. CONCLUSIONS 

We have performed HNC calculations of the free 
energy, internal energy, and pressure for various ionic 
mixtures with different fractional abundances of the ion 
species in a broad range of Z and T values. We have 
constructed an analytic approximation to the deviation 
from the LMR, which recovers the Debye-Huckel formula 
for multicomponent plasmas at T <C 1 and the LMR 
at T > 1, and which describes our calculations at any 
r values, as well as the HNC and MC results for the 
internal energy of plasma mixtures, available in the 
literature, with an accuracy better than 2%. 



The difference between the HNC results and formula 
(HHJ) for / ex lies within 0.013 and within 1%. The internal 
energy calculated using the analytic derivative of the fit 
([11]) deviates from the HNC results by not more than 
0.017 and not more than 1.5%. These maximal deviations 
are attained for the extremely asymmetric mixtures with 
Z 2 /Z 1 = 16. 

The analytic formula (|T6|) has also been compared to 
available HNC results for binary ionic mixtures in the 
polarizable electron background [J] and found to be sat- 
isfactory within the accuracy of the latter results (an ex- 
ample is shown in Fig. [J). 
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